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We present a second-order Godunov algorithm for multidimensional, ideal MHD. 
Our algorithm is based on the unsplit formulation of Colella (J. Comput. Phys. 
vol. 87, 1990), with all of the primary dependent variables centered at the same 
location. To properly represent the divergence-free condition of the magnetic fields, 
we apply a discrete projection to the intermediate values of the field at cell faces, and 
apply a filter to the primary dependent variables at the end of each time step. We 
test the method against a suite of linear and nonlinear tests to ascertain accuracy 
and stability of the scheme under a variety of conditions. The test suite includes 
rotated planar linear waves, MHD shock tube problems, low-beta flux tubes, and a 
magnetized rotor problem. For all of these cases, we observe that the algorithm is 
second-order accurate for smooth solutions, converges to the correct weak solution 
for problems involving shocks, and exhibits no evidence of instability or loss of 
accuracy due to the possible presence of non-solenoidal fields. 
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1 Introduction 



In this paper we present a new Godunov method for the equations of mul- 
tidimensional ideal magnetohydrodynamics (MHD). We give results from an 
implementation of the unsplit, second-order method of Colella [10] for these 
equations. The base scheme solves the ideal MHD equations using a second- 
order predictor-corrector formalism. To the base scheme we add three algo- 
rithmic components, whose effects upon accuracy and stability are measured. 
The first component is a MAC projection [6,19] step, which uses a Poisson 
solver to ensure that the cell-edge centered fields used to calculate fluxes are 
divergence-free to solver tolerance. This component, though commonplace in 
the context of incompressible Navier-Stokes simulation, is new to the MHD 
community. The second component is an approximate projection [2] that uses 
another solution of the Poisson equation to ensure that the cell-centered field 
is divergence-free to second-order. The last component is a filter [24] that also 
acts to suppress monopole sources in the cell-centered field. 

The section that follows covers recent work and some of the schemes used for 
ideal MHD simulation. It also introduces methods for enforcing the divergence- 
free constraint, and introduces some issues surrounding multidimensional MHD 
Section 3 introduces our basic algorithm and the extensions we have imple- 
mented. A suite of linear and nonlinear test problems will be used to deter- 
mine which of our algorithmic extensions are best suited to each problem type. 
These tests and results are covered in Section 4. The overall purpose is to find 
one combination of these extensions that is well suited to all of the problems 
considered. This will be done through comparisons to published results and 
in some cases to an eight-wave MHD algorithm we have implemented. We 
find that a projection step is indeed required for accuracy and stability of the 
schemes. For the eight-wave scheme in particular, use of the MAC projection 
is essential to obtain correct MHD shock jumps. 



2 Background 

The study of numerical algorithms for magnetohydrodynamics simulations re- 
mains an active one, with no one method having become the standard. Two 
generic algorithms are the most widely used at present: the Method of Charac- 
teristics/Constrained Transport (MOC/CT; common in the astrophysics com- 
munity) [15,33] and shock-capturing (Godunov) methods [5,9,11,12,18,30,36,39] 
Each has distinct benefits and drawbacks. Codes implementing the MOC/CT 
algorithm are relatively simple in design, and satisfy the divergence-free con- 
straint to machine precision. However, the method of characteristics used by 
the ZEUS scheme, as outlined in [33], is by construction second-order on Alfven 
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and advective waves, but does not address the two compressive waves of ideal 
MHD. Moreover, Falle [16] found that ZEUS exhibits spurious rarefaction 
shocks in certain 1-D MHD shock tubes for a non-isothermal equation of 
state. Codes implementing the shock- capturing algorithm on the other hand, 
while more complex, give highly accurate results even for strong shocks. They 
suffer from the drawback that the divergence-free constraint is only satisfied 
to truncation error, which can be large in the region of large jumps. In order to 
treat this difficulty, a variety of techniques have come into use. One such is the 
hybrid CT/shock-capturing scheme [5,13,23,31,36], for which the constraint is 
satisfied by design like in the MOC/CT case. The cost to the accuracy of the 
underlying shock-capturing scheme is unclear. Another approach, originally 
due to Brackbill and Barnes [8] and implemented by workers such as Ryu et 
al [30], uses a divergence cleaning step on the cell-centered fields to enforce 
the constraint. In a third approach, Powell and co-workers [18,27] use a eight- 
wave reformulation of the ideal MHD equations originally due to Godunov 
[17]. Toth [36] (hereafter TOO) implements all three types of schemes, among 
others, using them as the basis for a comparison on a variety of 1-D and 2-D 
tests. More recently, Dedner et al [14] compare several hyperbolic schemes 
with additional waves and divergence-damping terms on the TOO tests. 

The 2-D tests of TOO serve to underscore the importance of using multidi- 
mensional problems in evaluating different algorithms, since it was mainly in 
this context that differences between them became apparent. This is to be 
expected, since errors due to non-solenoidal fields will generally only show up 
for problems in two or more dimensions. Many shock- capturing MHD schemes 
use an operator-split (or dimensionally- split) formalism to treat multidimen- 
sional ideal MHD. This means that, for each spatial dimension of the scheme, 
the one dimensional MHD equations are applied once. Unsplit schemes, which 
instead use the full multidimensional version of the equations, have been im- 
plemented and shown to give results equivalent to those of split methods for 
hydrodynamical problems [10]. Unsplit shock-capturing schemes for multidi- 
mensional ideal MHD are relatively new, however. One of our main goals is 
to assess the efficacy of different approaches for enforcing the divergence-free 
constraint in one such unsplit scheme. 



2.1 The Divergence-free Constraint 

Some means must be employed to ensure that the field satisfies the divergence- 
free constraint, since this is only guaranteed to within truncation error in 
shock-capturing schemes. A dramatic example of even small errors in the 
divergence-free condition leading to instability is given in Section 4.1. There, 
small amplitude, non-propagating waves become unstable unless non-solenoidal 
fields are smoothed. This phenomenon can be explained in terms of a modified 
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equation analysis, discussed in Appendix A. The possibility of incorrect field 
topologies, incorrect dynamics, and numerical instability motivate efforts to 
formulate and understand different means of supressing non-solenoidal errors 
in the magnetic field. 

The eight-wave MHD algorithm, as implemented by Powell et al [18,27] and 
others, addresses the problem by adding additional terms corresponding to 
monopoles to the ideal MHD equations. The resultant equations are sym- 
metrizable, so that they are Galilean invariant and transport V • B [17]. The 
additional terms show up in two ways for shock-capturing schemes. Since they 
modify the 1-D MHD equations used for characteristic tracing to include an 
additional eighth wave that travels at the flow speed, monopoles will be ad- 
vected along with the flow. Such monopoles could be carried out of the domain, 
or they might build up at a stagnation point. Secondly, the additional terms 
appear as source terms, making the system non-conservative if the divergence- 
free condition isn't already satisfied. 

Dedner et al [14] test a scheme that extends the eight-wave concept to damp 
and advect monopole sources, even at stagnation points in the flow. This 
is done through the magnetic analog of an artificial compressibility term, an 
approach that appeared earlier in the context of the Maxwell equations [24]. To 
the extent that the computed auxiliary field remains continuous, the scheme 
will remain divergence free. Any monopole sources will be advected at the 
fastest speed allowed under the Courant condition, and damped as they are 
advected. This method is very useful on unstructured grids, where solving the 
Poisson equation in order to project out the solenoidal component of the field 
is difficult. 

Divergence cleaning, or Hodge projection, in shock-capturing schemes can ad- 
dress the problem of non-solenoidal fields in two ways. The most widely dis- 
cussed [3,8,30,36] involves projection of the cell-centered field onto the space of 
divergence-free fields. Projecting in this manner with a centered difference ap- 
proximation to the divergence is consistent with the underlying cell-centered 
scheme. One is left with fields at the advanced time which are divergence- 
free to machine accuracy. Such a projection has been found to give correct 
field topologies in shock tube problems [30]. On the other hand, choosing to 
eliminate the divergence as measured in one metric does not guarantee that 
unphysical effects are not entering into the dynamics. An illustrative example 
from incompressible flow [22] where an analogous constraint on the fluid veloc- 
ity, V • v — 0, occurs shows that checkerboard modes in the velocity can cause 
instabilities when using a centered difference to approximate the divergence. 
Such modes must be damped by a suitably chosen filter in order to regain 
stability. 

Another option for cleaning of non-solenoidal fields is to project the fields at 
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cell-edges, which are then used to calculate fluxes. In a MAC discretization 
[6,19], the divergence-free constraint is enforced to within solver tolerance, 
which can be set as low as machine precision. We use a multigrid solver to 
solve the associated Poisson equation. The MAC projection has the advantage 
that it does not affect the conservation properties of the scheme. 



2.2 Multidimensional MHD 



The MHD equations in two or more dimensions are decidedly more complex 
to solve than in one dimension. In MHD simulations with variations along the 
x-axis alone, there is no change in the field along the x-axis. The divergence- 
free constraint is therefore trivially satisfied. Obtaining the solution to the 
Riemann problem upon which Godunov methods are based is also relatively 
straightforward in 1-D. 

In multiple dimensions, we are solving the full equations of ideal MHD, which 
in conservation form are 



d t { P u) + v ■ 



( B 2 \ 1 ^ 



d t p + V-{pu)=0 




uB - Bu 



d t ( P E) + V 



pE + P + Ib>)u-±-(u.B)B 




0. 



(1) 

(2) 
(3) 
(4) 



— * 

subject to the constraint V ■ B = 0. Here p is the mass density, pu the mo- 
mentum density, B the magnetic field, and pE = \p\u\ 2 + g^|-B| 2 + t~[P the 
total energy density. The d t notation denotes derivatives with respect to time. 



3 Equations and Algorithms 



In general for hyperbolic conservation laws, the conserved variables U evolve 
according to dtU + V • F(U) = 0. Our code uses a second-order Godunov-type 
method for hyperbolic conservation laws [7]. Such volume-average schemes 
follow the flux of conserved quantities such as momentum into and out of 
each cell comprising the computational domain. Quantities are stored at cell- 
centers, their value at this point being the volume-average over the entire cell. 
During the course of a timestep, the flux of the conserved quantities at each 
edge of each cell is computed. Differencing these fluxes gives the update of the 
conserved quantity to the next time (see Figure 1). 
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Fig. 1. Single cell of the computational domain, showing interpolated states W i • ± x 
and W i j^ y , and cell-centered state Wij. The (Riemann) states at cell edges would 

n+i 

lie between the interpolated states. So W ± , for instance, would lie between 

<t and 

More specifically, our code is based on the unsplit MHD algorithm of Colella 
[10]. It ostensibly tracks the evolution of the conserved density p, three com- 

— * 

ponents of momentum pu, three components of magnetic field B, and total 
energy density pE, in two spatial dimensions. These variables would then be 
evolved according to the ideal MHD equations (1-4). However, for simplicity 
and accuracy a set of primitive variables W consisting of the density, velocity, 
field, and pressure are evolved in time: 

W = [p u v w B x B y B z P\. (5) 



We will switch freely between U and W here depending on which is most 
convenient; the indices will indicate the centering of the variables. Thus i,j 
will denote cell-centered states, ± forward- and backward-interpolations 
to cell-edges, and i + |, j and + ^ (Riemann) states at the edges. 

The fundamental aspects of our scheme are as follows. Given the state Wfy at 
time n and spatial coordinate (xi,yj), we simultaneously interpolate in space 

and extrapolate in time to obtain the states W i; j^ x on the two x-boundaries 

of each cell The same is done for the states at the y-boundaries, W { j ^ y . 
All this is done in the normal and transverse predictor steps of the algorithm. 

~ n+i 

Next, the Riemann problem is solved using the W i j £ x states. These Riemann 
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states are used to calculate fluxes at each edge of a cell in the corrector step. 
The fluxes are then differenced to give the update to the next timestep, W™^ 1 . 
We go into more detail on these steps of our scheme in Sections 3.1 through 
3.4. 

The additional parts of our MHD scheme consist of additional terms in the 
normal predictor step to ensure correct multidimensional behavior, and several 
steps that address the divergence-free constraint. These are also discussed in 
Sections 3.1 through 3.4. 



3. 1 Normal Predictor 



The predictor computes W™*± x and W™*± y using a Taylor series expansion: 
<1 = W Z ± %0,W - ^A*8,W - ~d»F>. (6) 



Here the matrix A is related to the flux F x of of the conserved variable U 

by A = (duW) (dwF x ). We give formulae for W^± jX , those for W™*£ y being 
similar. The first three terms on the right-hand side are computed in the nor- 

mal predictor, and we label this intermediate result W^- ± x . We now separate 
out the evolution of the normal field B n = B x through the following notation: 
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w = 




, A = 




B x 








(7) 



The matrix A corresponds to the usual 1-D MHD equations, with its seven 
characteristics: forward- and backward-propagating fast, slow, and Alfven, 
plus the advective wave. Second-order accuracy is achieved in part through 
the use of characteristic analysis to calculate second-order accurate derivatives 
d x W in the spatial interpolation. This characteristic interpolation is based on 
calculation of the eigenvalues and left- and right-eigenvectors Ik and of 
the matrix A, giving the following expression for the interpolation of the W 
variables to cell-edges: 



h3 



E 



At , 



OLkTk 



Oik 



Min(a°, a + ,a~) if a + a~ > 
otherwise 



(8) 
(9) 
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a- = 2/ fc • (Wft - W7L,,.) 



(10) 

(11) 
(12) 



The ttfc represent the strength of the k th wave in the interpolant. The sum 
over \ k < would correspond to backward-propagating waves used in the 

interpolation to W^^ x , and similarly for X k > and H^ j|+ X . 



A full accounting for all x-derivative terms in the 2-D MHD equations shows 
that as is given by 



a B = 



0, 



B T B„ 



B 7 



Airp Awp A7vp 



u-B 

47T 



(13) 



These terms are essential to the second-order accuracy of the scheme, in par- 
ticular on multidimensional problems such as waves not propagating along the 
coordinate axes. They are incorporated into our algorithm through a simple 
finite differencing of the normal derivative of the normal field. The terms are 
added to those already present due to the characteristics-based interpolation: 



rv l,J,±,X 



W^l x -^a B (V x B:) h3 . 



(14) 



The (V x B%) itj = ((B x )? +1 j - (B x )^_ 1J )/(2Ax) correction term is the centered- 
difference approximation to d x B x . Note that this approximation to the deriva- 
tive is not limited. The need for such correction terms in an unsplit scheme for 
multidimensional MHD was not addressed in C90. It was first noted by Stone 
[34] during an examination of the accuracy of an unsplit Godunov scheme on 
advected flux rings. 



3.2 Transverse Predictor 



The last term in the evolution equation (6) is included via the transverse 
predictor. The basic idea is to approximate transverse derivatives (in this 
case in the y-direction) using a 1-D Godunov method. We take the states 

calculated in the normal predictor, W^^y, and first use them to solve the 
Riemann problem at each y-boundary in the domain. The resulting Riemann 

states U \ are subsequently used to calculate the fluxes needed for the last 

i,j+2 

term in equation 6. In more formal terms, 
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u 




i,3+2 




= F y [U 




) 



) 



(16) 



(15) 



Here 1Z(-, •) denotes the Riemann problem solution using the two states on 
either side of an edge as input; see Section 3.4. It is then a straightforward 
matter to use these fluxes to calculate a finite-difference approximation to 
d y F y , and thereby complete the calculation of the edge-centered states. 

3.3 Corrector 

The corrector first calculates fluxes at all cell-edges using another Riemann 
problem solve. At this stage, we have the first-order accurate approximation 

to the interpolated states j ± . from the predictor in hand. The Riemann 
solver takes these states and returns a single state for each cell-edge, U x 2 



The formula for U \ is similar. See Section 3.4 for more details on the 

i,j+2 

Riemann problem solution. 

The states at cell-edges are not guaranteed to be divergence-free. We modify 
the C90 algorithm to enforce the divergence-free condition for the Riemann 
problem states. These Riemann states have a non-solenoidal component that 
we treat through a MAC projection, earlier used in the context of incompress- 

ible fluid computations [6]. The edge-centered fields (B*) 1 and (B*) 2 1 

are used to calculate a cell-centered monopole charge density qu = V • B* , 
and a Poisson solver is in turn used to find the scalar field implied by this 
monopole charge density distribution. The scalar field satisfies the following 
relations, in which T>^ correspond to the forward- and backward-difference 
approximations to the derivative and similarly for 



and U \ . For instance, 

M+2 




(17) 
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(18) 
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The correction to the field is calculated from 6 as follows: 



(19) 
(20) 









< B »>, + i, 




i 

2 








1 

2 











(21) 

(22) 

(23) 
(24) 



With this correction to the magnetic field of the Riemann states, the L 1 norm 
of the MAC monopole density is reduced from its initial value by a user- 
settable multiplicative factor, in our case 10~ 12 . 

The algorithm now proceeds to calculate the fluxes associated with the Rie- 
mann states. These fluxes are then differenced to give the update to the next 
time mt 1 : 



1 / 1 



F \ 2 =F X [U \ 2 (25) 

*+2 J V J +2'V 



= UTj - At V-F X 'f 2 - At V-F V ; n ^ (26) 

After the update to t n+1 , we are left with a cell-centered field B*f +1 that 
is no longer divergence-free by a centered-difference divergence metric. To 
what extent, if any, this is a problem depends on the physical problem being 
considered, and will be addressed later. 

For those cases where a reduction of the divergence is required, two algorithmic 
extensions have been implemented. The first follows from noting that it is 
desirable to have a diffusive term of the form (see [24,14]): 

= »V*(V ■ B) (27) 



act on the divergence of B. This may be rewritten to eliminate a spatial 
derivative by pulling out a divergence operator, giving 

— = r?V(V ■ B) (28) 
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A simple, single-step filter may be derived [26] as a finite-difference approxi- 
mation of equation 28: 



B x := B x + V At (VtV-B x + V° x V° y B y ) (29) 
B y := B y + V At (V° x V° y B x + VfV-By) (30) 

The advantage of this choice of discretization of Equation 28 lies in its effec- 
tive damping of checkerboard modes; see [28] for analysis and comparison of 
different filtering schemes. 

In order to choose a value for rj we use Fourier stability analysis, giving a 
stability condition for the scheme in equations 29 and 30 of At < . Since 
At is set by the Courant condition, we are able to derive a condition on rj, 
giving the maximum amount of diffusion of the monopole sources possible, 
given our timestep and grid spacing: 



with C < |. A stronger condition, C < |, will always damp monopole modes. 
This formulation will both decrease the cell-centered divergence and damp 
checkerboard modes. The filtering is always used when an approximate pro- 
jection is performed. 

The numerical effect of the filter can be modified through the parameter C . 
We found, through the nonlinear tests outlined in Section 4, that values in the 
range 10~ 2 — 10 _1 worked best. We note here that we chose to apply the filter 
to the conserved variables U n+1 , so that the changed magnetic field causes 
no change in the total energy. Any addition (subtraction) of magnetic energy 
shows up as a decrease (increase) in the internal energy. 

The second means for treating non-solenoidal cell-centered fields is an approx- 
imate projection, in which the cell-centered divergence is used to calculate a 
monopole charge density qu = 'D x , B* + V y ) B*. Note that the multigrid Poisson 
solver uses the finite-difference operator V + V~, not V°V°, for the Laplacian 
V 2 . The resulting solution is not discretely divergence-free, as it would be if we 
used the operator V°V° [2], but instead second-order accurate. However, the 
extension of exact solvers to adaptive meshes is extremely complicated [20], 
whereas the approximate projection used here is straightforward to extend 
to AMR [25]. We note moreover that use of a centered-difference measure 
of divergence leaves unchanged unphysical checkerboard modes [22,26]. An 
extensive analysis of approximate projections is given in [1]. 

The Poisson solve yields the scalar field satisfying equation 19. We difference 
as follows to give the corrected field: 
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B n + 1 = BT+ 1 _ V 0^ . ( g 2) 

B n y +1 = B* y n+1 - (33) 

We find in our tests that an approximate projection does not have a significant 
impact on the quality of our results when used in conjunction with a MAC 
projection. Given the significant computational expense of using both, we 
advise against such a scheme. 



3.4 Riemann Solver 



We solve the Riemann problem for ideal M HD using a linearized solver [35] . It 
employs characteristic analysis, like that of the normal predictor of Section 3.1, 
to solve for the states at cell edges. In this case, eigenvalues \ k and eigenvectors 
l k and r k are calculated using the arithmetic average of the states to the left 
and right of the edge. (Call this the base state.) For example, the solution to 
the Riemann problem at x. 1, W RR , is built from the appropriately directed 



i+2 



jumps across waves as follows: 



W RP = { 



\ (W R - Ek,x k >o ®kr k + W L + Ek,x k <o ®kr k ) , if |w| < 

W L + E fc ,A fe <o ®kr k , if u < -e v (34) 
W R - Efc,A fc >o ®kr k , if u > e v 



Here \u\ is the advective velocity (along e x in this case) of the base state. The 
constant e v is chosen to be 10 -14 times the largest characteristic speed X k . 
Note the averaged solution used for small advective speeds. This ensures that 
advective velocities on the order of machine precision, whose sign is random, 
do not cause asymmetries in the Riemann problem solution. Finally, note that 
in order to maintain consistency, the longitudinal component of the magnetic 
field in the solution (in this case B RP ) is assigned the value of the longitudinal 
magnetic field of the base state, \{B% + B R ). 

For strongly nonlinear problems, it sometimes happens that the CFL condition 
is not sufficient to keep a scheme stable. Ordinarily, we use a simple minimum 
of timesteps implied by the one-dimensional CFL condition. If accelerations 
are large, velocities can grow such that in one timestep the pressure or den- 
sity becomes negative. In order to dynamically adjust to such situations, we 
implement a scheme that checks for negative pressures or densities in the cell- 
centered states, W™^ 1 , after the corrector step. If one is encountered, that 
timestep is restarted with the CFL number lowered by a factor of two from its 
nominal value, down to a minimum of one-eighth of the initial CFL number. 
Once the code has proceeded for several timesteps without again encountering 
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negative values, the CFL number is raised by a factor of two, eventually reach- 
ing its initial value. In practice, negative values are encountered only rarely, so 
that the average CFL number over the course of a run is close to the nominal 
value. Note that even if the CFL ramps quickly back up to the nominal value, 
the timestep is never allowed to increase by more than 10% in one iteration. 
In this way, the effect of an increase in CFL number is in fact spread over 
several iterations. Note that the average CFL number remains above 0.4 in all 
tests of Section 4. 



4 Numerical Tests 



In this section, we compare the behavior of the code on a variety of linear 
and nonlinear problems. Both seven- and eight-wave MHD codes were used. 
The code implementing the eight-wave MHD algorithm also uses a predictor- 
corrector formalism. The characteristic analysis performed in the predictor 
steps uses an eighth wave carrying changes in the normal field at the advection 
velocity, in addition to the seven waves of ideal MHD. As a result, in equation 
7 of Section 3.1, the (2,2) entry of A is equal to the advective speed u and 
not zero. The terms needed for second-order accuracy in ideal MHD (equation 
13) are accounted for in this case, so that as = 0. The eight-wave algorithm 
also adds a source term, 



S = 



Bx By_ Ih u-B 

47r' 47r' 47r' ' 47T 



(V-B), 



(35) 



to the right-hand side of the MHD equations 1-4. The source term is calculated 
in the transverse predictor and corrector steps. It is added to the updated 
states along with the differenced fluxes. We note, however, that performing a 
MAC projection guarantees that V ■ B = for the edge-centered states, and as 
a result the source term calculated in the corrector is numerically zero in this 
case. This implies that only with a MAC projection is the scheme conservative 
and satisfies the jump relations. The source term in the transverse predictor, 
while generally small, is not numerically zero. 

For both MHD implementations, different variations on the base algorithm 
were tested. In what follows, codes with conservative filtering are labeled with 
'CF'. An 'AP' denotes codes with an approximate projection; it is always ac- 
companied by conservative filtering, in order to suppress checkerboard modes. 
Those codes utilizing a MAC projection are labeled with 'MAC. So, for ex- 
ample, a MAC projection code with filter is labeled 'MAC+CF'. 
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4-1 Simple Linearized MHD Waves 



We tested the performance of the code on all four varieties (advective, fast, 
slow, and Alfven) of linearized MHD waves. These tests comprised waves prop- 
agating along x- and y-axes (wavenumbers n = (1,0) and n = (0, 1)), along 
with waves at slopes of 1:1 (n = (1, 1)) and 2:1 (n = (2, 1)). The simulation 
domain had length L — 1 in both dimensions, and the boundaries in all cases 
were periodic. The Alfven waves used p — 1, u — 0, B — B b = y/4arb with 
unit vector b — (^, ^j), and P — 1. The perturbation 5W is 



5W 







-c A 



Bo 




5 pcrt sin (k ■ xj , with k = 2nn, 



(36) 



where n is the aforementioned vector of integers chosen so as to be consistent 
with the periodic boundaries. The Alfven speed ca = B / ^Airpo = 1. 

Fast and slow wave expressions are somewhat more complicated. In this case, 
b lies at 45 degrees from the unit wavevector k, all other aspects of the un- 
perturbed state remaining unchanged from the Alfven case. The perturbation 
is 



8W = 



Po 



CF/S 



a 2 h x -^2c 2 p/s b 



C F/S 





5 

A 

C 2 -a 2 



-n,, 



-n. 



p a 



5 pert sin (k ■ xj (37) 
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a 2 = ^ (38) 
Po 

c f/s = 2 ( fl2 + 4^ + a 4 ^ corresponding to fast/slow speeds (39) 

71 

n = (n x , fly) = — — (40) 

77/ 

All tests gave results consistent with second order accuracy, although waves 
at 2:1 showed slightly smaller (> 1.8, versus 2.0) rates of convergence in some 
components. The modifications suggested by Stone were essential in obtain- 
ing second-order convergence for the 1:1 and 2:1 slope tests, increasing the 
convergence rates from first- to second-order for waves propagating along the 
diagonal. An example of the convergence rates for small amplitude fast waves 
is given in Table 1. 



Component 


Without correction 


With correction 


P 


0.977 


2.03 


V x 


1.40 


2.03 


Vy 


1.15 


2.02 


B x 


1.60 


2.01 


By 


1.57 


2.02 


P 


0.977 


2.03 



Table 1 

Convergence rates by component for fast waves (<5 per t = 10~ 5 ) propagating at 45 
degrees, for seven- wave MHD code without and with corrections suggested by Stone 
[34]. 

A more stringent test of the code is to advect the linearized waves so that their 
profile remains stationary. In these tests, the background state has a non-zero 
velocity equal to the wavespeed: 



Alfven: u = —caU (41) 
Fast/Slow: u = —CF/sn (42) 

An analysis of the eigenstructure of MHD shows (Appendix A) that such waves 
should cause trouble for the seven-wave MHD codes that do not suppress non- 
solenoidal fields. Errors generated in this case are not advected away, causing 
difficulty in the case that they are not diffused or otherwise dealt with. We find 
that the seven-wave MHD algorithm without the application of any projection 
or filter is unstable for low amplitude (S pert = 10~ 4 ) fast waves with n at 2:1 
slope. The instability starts as a high frequency oscillation in the field parallel 
to n. It grows with time and spreads to the other components, causing a low 
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order of convergence even at early times. Once the oscillations reach a certain 
level, the solution becomes unstable. The oscillations are almost completely 
absent, and the code stable, using a MAC projection. Adding filtering of the 
field further reduces the errors. A code including filtering but no MAC pro- 
jection, and one with approximate projection and filtering both stabilize the 
scheme. We conclude that either projection or filtering is essential to stability 
for this class of problems. 



4-2 Decay of Linearized MHD Waves 



The equations of ideal MHD neglect the effects of viscosity and electrical 
resistivity. Numerical dissipation, however, can affect the solution in ways that 
mirror these physical effects. Following Ryu, Jones, and Frank [30] (hereafter 
RJF95), we measure the decay of Alfven, fast, and slow waves, and use the 
implied physical resistivity as a measure of the numerical resistivity of our 
ideal MHD scheme. 

We use exactly the same set-up as RJF95, with Alfven, fast, and slow waves 
propagating at a 1:1 slope and wavelength of y/2 times the length of one side of 
the computational domain. The decay of these waves was measured by fitting 
a decaying exponential to a measure of the wave strength, 

£ = E'*-Wj. ( 43 ) 



where (5W)ij = Wij — Wq is the perturbation in the primitive variables. 
is the left-eigenvector associated with the mode in question, evaluated at Wo, 
the unperturbed state. Figure 2 shows the dependence of the decay rate on 
resolution for runs with number of cells per dimension N =16, 32, 64, 128, 256, 
and 512. Comparison with Figure 4 of RJF95 reveals a much smaller decay 
rate for our unsplit method. Furthermore, we find that the decay rate varies 
according to the power law T oc A^ -3 , not N~ 2 . The former is consistent with 
second-order accuracy, since the truncation errors have the following form: 

- = ^(Ax) 2 ^ + ^(As) 3 ^ + 0{Ax% (44) 



We see here that the second-order error term is proportional to ^Jf, which 
is dispersive in nature. The third-order term is the first error term that is 
dissipative. As a result, we expect the dissipation to decrease as Ax 3 , or A^~ 3 . 
We have no good explanation for the 0(Ax 2 ) decay law observed in RJF95, 
aside from the possible effects of artificial viscosity. 
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Fig. 2. Normalized damping rate versus resolution for the decay of Alfven, fast, and 
slow modes propagating at a 1:1 slope. The calculations were done at resolutions of 
16 x 16, 32 x 32, etc. up to 512 x 512. 

4.3 MED Shock Tube 



The second test is the MHD shock tube problem from Ryu, Jones, and Frank 
[30] . The solution consists of two fast shocks, one slow shock, one slow rarefac- 
tion, and a contact discontinuity. We have run the problem in two orientations: 
(1) with the shock velocity aligned with the x-axis of the computational do- 
main (referred to below as 1-D), and (2) with this velocity inclined at a 2:1 
slope. The latter configuration follows the 2-D shock tube test case from TOO. 
It was chosen in order to test the multidimensional behavior of the code, mean- 
while ensuring there were no serendipitous cancellations of errors, as might be 
the case in a 45° inclined case. In runs of a 1-D tube at i? 5 i 2 (ie. 512 cells per 
linear dimension) on a domain of size L — 1 to a time of t — 0.08, we are able 
to reproduce the results given in Table VI of Dai & Woodward [11] (hereafter 
DW94) with errors of < 0.12%. These results are independent of the type of 
filter used and whether an approximate projection was performed. Both the 
seven- and eight-wave codes give exactly the same result. When combined 
with an observed first-order convergence rate, they give us great confidence in 
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the performance of both codes on 1-D shock tube problems. 

The inclined version of this shock tube problem was run on grids of size 2NxN, 
with N = 64, 128, and 256. Unlike in the aligned shock tube, in this case 
divergence cleaning is required. Since the natural boundary conditions for 
the problem would be tedious to implement in our multigrid solver, we chose 
instead to embed the physical domain within a larger one, four times its size. 
The data in the region outside of the original one is filled using the natural 
symmetry of the problem. The Poisson equation is then solved on this enlarged 
domain, using homogeneous Neumann boundary conditions. We have verified 
that the field in the original domain is clean to within the tolerances used 
previously. 

We compared the results from the inclined shock tube with coarsened versions 
of a -R4096 1-D shock tube results by first taking cuts of the data that included 
all cells lying along a line at 2:1 slope (a = 26.57°). Then, all velocities and 
fields in the cut were rotated by —a. A check using the initial conditions 
showed that such a cut of the inclined initial conditions matched perfectly 
with the 1-D shock tube ICs outside of the jump region. (In the jump region, 
the inclined ICs were slightly different due to the volume averaging performed 
in producing them. This leads the initial jump to be spread over two zones 
instead of one.) Note that we chose the grid spacing for the inclined runs to be 
a factor of \fh smaller than the 1-D grid spacing, so that the shock covers the 
same distance in physical space in a given time. The domain was therefore of 
length L = 1 in the direction of shock propagation, and the final time t = 0.08, 
in both the aligned and inclined shock runs. 

Figure 3 shows results for the inclined shock tube overlaid on the (coarsened) 
aligned result. They compare a i?256 1-D run with a N = 256 2-D run that has 
equivalent effective grid spacing. The ideal MHD codes shown give fractional 
errors with an L 1 norm of 3% or less, averaged over all components. Errors in 
the normal component of the field are relatively large, however. As illustrated 
here and by Figure 11 in TOO, shock-capturing codes generally produce un- 
physical variation in the normal component of the magnetic field inside the 
shock transition layer in non-grid-aligned shock tube problems such as this. 
More generally, such codes exhibit non-monotonic behavior inside the shock 
structure, so that Riemann invariants are not exactly preserved there [38]. 
Importantly, despite errors in B n inside the jump region, the jump relations 
are still satisfied outside of it. This is clearly illustrated in Figure 4, where we 
see the normal field for seven- and eight-wave schemes with filtering. Note in 
particular that while the eight-wave scheme converges to an answer, it does 
not give the correct shock structure. 

Figures 6 and 7 show plots of the L 1 norm of the error versus resolution for 
the 2-D shock tube runs. We first note that neither the seven- nor eight-wave 
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Fig. 3. MHD shock tube comparison by component. The 1-D shock tube (line) 
is used as a basis for comparison with 2-D shock tube runs with same effective 
resolution. The jumps in normal field B n are expected, and the size of the jumps is 
similar to that published in Toth 2000. 
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Fig. 4. Normal field results for MHD shock tube problem. Resolutions from Rq4 to 
i?5i2 are shown. Note the incorrect shock structure in the eight-wave result. 

base code converges at first order. This is most evident in the errors in the 
transverse velocity and transverse field. The same plots indicate filtering alone 
is not sufficient to produce reasonable convergence rates. A MAC projection 
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step added to either scheme produces much improved results (label 'MAC 
in the figures), and the further addition of a filtering step gives first-order 
convergence in all components (label 'M AC+CF'). The use of an approximate 
projection (labels 'CF+AP' and 'MAC+CF+AP') produces decreased error 
in the normal field, indicating a more accurate solution at given resolution. 
However, no improvement in the convergence rates or errors in the other five 
components is observed. When averaged over all components, the decrease in 
error observed with an approximate projection is still noticeable, though less 
so than indicated by considering the normal field alone. 

While the results for seven- and eight-wave codes with a MAC projection step 
are very similar, this picture changes somewhat when the MAC projection is 
replaced with an approximate projection (label 'CF+AP'). The seven-wave 
code gives identical results whether or not a MAC projection is used with the 
approximate projection and filter. The eight-wave code still exhibits smaller 
errors when an approximate projection is utilized. However, convergence of 
the errors for the eight-wave code with approximate projection alone stalls 
at higher resolution. The reason for this behavior becomes more clear upon 
closer inspection of the normal field, plotted in Figure 5. The normal field 
converges to an answer that is 0.1% above the expected value. This error, 
while masked by the larger errors at the jumps at lower resolution, becomes 
significant at higher resolution and consequently hampers convergence. This 
error is observed to grow steadily over time. 

4-4 Magnetized Flux Tube 

This problem involves a high-field, low gas pressure region bounded on both 
sides by a high-gas pressure, zero-field region. The physical domain is of length 
L x = L y = 1 on both sides. The base state has the entire 2-D domain in pres- 
sure balance. The boundaries between magnetized and unmagnetized regions 
are discontinuous and lie along x = ±0.2. In both regions, p — 1, u — and 
B x = B z = initially. In the magnetized region, x E [—0.2,0.2], B y = ^/80n 
and P — 1. Outside of this, B y = and P = 11. To the background state 
is added a sinusoidal perturbation upon the x-velocity whose amplitude is 
£ P ert = 0.01 times the Alfven speed, and covers the entire domain: 

5u = 5 pcrt CA sin (2ny) and SB = 5p = 5P = Sv = 5w = (45) 

The Alfven speed in this case is ca = B> / \JAnp = \^20. The strong discon- 
tinuity in the field is expected to cause problems for algorithms that do not 
suppress non-solenoidal fields. In particular, at the stagnation points where 
the perturbation velocity is zero, truncation errors leading to non-solenoidal 
fields can build up and cause numerical schemes to go unstable. This is in fact 
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Fig. 5. Normal component of the field for eight- wave codes with (bottom) and 
without (top) a MAC projection step. In these results for the inclined shock tube, 
a cut along the shock propagation direction is shown. The eight-wave code with the 
MAC projection converges to the correct answer, B norrna i = 5. Without the MAC 
projection, the code converges to an incorrect value for the normal field, about 0.1% 
higher than expected. 
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Fig. 6. L 1 norm of error versus resolution (number of cells per dimension) for 2-D 
MHD shock tube run with variants of the seven-wave code. The dashed line is a 
fiducial showing first-order fall-off of errors. Note that some manner of projection 
step is required for first-order convergence. 
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what we find for base schemes run without a projection step. Non-solenoidal 
fields build up at these stagnation points, causing the production of spurious 
velocities and triggering numerical instabilities. 

Figure 8 shows the initial conditions for the flux tube problem. The pertur- 
bation plucks the field. After an initial transient at start-up, the perturbation 
develops into a standing wave in the magnetized region. The unmagnetized 
region sloshes back and forth from one side of the flux tube to the other, owing 
to the periodic boundary conditions. We expect the initial velocity perturba- 
tion to cause standing Alfvenic waves in the magnetized region and, because 
of the slightly different oscillation frequencies of the two regions, compres- 
sive waves at the boundaries of the tube. The change in the internal energy 
tracks the compressive waves. The problem was run to a time t = 6.0, cor- 
responding to about 60 Alfven crossings of the short dimension of the tube. 
This was enough time for the problem to exhibit stable periodic behavior and 
subsequently evolve for many periods of this oscillation. 

In this case we tested both seven- and eight-wave codes, but only with either 
a MAC projection and filter, or an approximate projection and filter. This is 
justified by the sub-par performance of the other variants on the inclined shock 
tube problem. These four versions of the code performed almost identically 
on this test, as evidenced by Figure 9. There, we plot global quantities such 
as kinetic and magnetic energy and L 1 norm of the monopole density. Note 
that the cell-centered divergence of the field is smaller for the approximate 
projection, though the dynamics remain the same. 



4-5 Inclined Flux Tube 

A second version of the flux tube problem, in which it is inclined at an angle 
of 45 degrees with respect to the original, is better at differentiating between 
the algorithms. It constitutes a strong test of the robustness and stability of 
the codes. 

The size of the domain was L x = L y = \^2, and both magnetized and unmag- 
netized regions have the same physical extent as in the aligned case. In both 
regions, we have p — 1, u — and B z = initially. In the magnetized region, 
B x = — V407T, By = V407T, and P = 1. In the unmagnetized region, P = 11 
and B x = B y = 0. The perturbation is again applied to the entire domain, 
and has strength <5 pe rt = 0.01: 




(46) 
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Fig. 8. Grid-aligned (top) and inclined (bottom) flux tube initial conditions. Vectors 
indicate the perturbation velocity field; it is constant along lines parallel to the 
vectors. There are two regions: one magnetized, the other not. The entire domain is 
initially in total (magnetic plus thermal) pressure balance. The velocity perturbation 
at 1% of the Alfven speed is applied to the entire domain, causing the tube and 
surrounding unmagnetized medium to oscillate. 
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Fig. 9. Grid-aligned flux tube results at #512- We use the L 1 norm of V • B, plus 
total energies, to follow the evolution. Energies are normalized to their initial values. 
All quantities are plotted in units of the Alfven crossing time, t a = l/c a = 0.2236. 
Despite their different means of suppressing monopoles, the codes performed very 
similarly on this test. 
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5v = 5 pcrt - 7 =sin[2n— /r j 
SB = 5p = 5P = 5w = 




(47) 



(48) 



One effect of rotating the flux tube is the addition of a numerical perturbation, 
or gridding effect. The volume-averaging required for producing the initial 
conditions creates a region of intermediate pressure and field between the fully 
magnetized and unmagnetized regions. The total pressure in this intermediate 
region was kept the same as elsewhere in the domain. The physical extent, and 
presumably physical effects, of this transition region become smaller as the 
resolution is increased. The second effect of our rotation is that the physical 
domain was somewhat larger. As mentioned above, this was done in order to 
keep the physical width of the magnetized and unmagnetized regions the same 
as for the aligned flux tube. The result is a factor of \[2 mismatch between 
the grid spacing in the aligned and inclined flux tube runs. 

We expect many qualitative similarities between the aligned and inclined flux 
tube results. The characteristics of the start-up transients are somewhat dif- 
ferent, due in part to the gridding effects. The slightly different grid spacing 
also means that slightly higher resolution is required to obtain comparable 
results. Indeed, these tests were run at resolutions ranging by factors of two 
from i? 12 8 to -R1024, a factor of two higher than for the aligned case. 

We again chose to run only those codes that showed close to first-order con- 
vergence for the inclined shock tube on this problem: seven- and eight-wave 
codes with either MAC projection and filter, or approximate projection and 
conservative filter. There was very little difference between the two seven-wave 
codes. The difference between seven- and eight-wave codes was more notice- 
able. The eight-wave codes both experienced large increases in kinetic energy, 
indicating that the magnetized tube was no longer in equilibrium. The size of 
this increase was largest (~ 10 3 KE in i t i a i) for the eight-wave code with approx- 
imate projection, which also experienced substantial deviations from energy 
conservation. The size of these kinetic energy jumps decreased with resolution. 

Figure 11 shows a resolution study of the inclined flux tube, with global quan- 
tities plotted out to t — 1. (This is about 5 Alfven crossing times; a shorter 
time was chosen in order to show more detail. All codes that were stable to 
this time were also stable out to t — 6.) The plots include kinetic, internal, and 
magnetic energies, all scaled to their initial values. We also show the L 1 norm 
of cell-centered measure of the divergence, D° ■ B, along with the magnetic 
energy change due to the MAC projection and filtering steps. First, note that 
the internal and magnetic energies appear to be converging, while convergence 
of the kinetic energy is not so clear. This result can be attributed to the large 
difference in the size of these quantities; the initial kinetic energy is ~ 10 5 
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Fig. 10. Inclined flux tube results for one eight- wave ("8+MAC+CF") and two 
seven- wave ("7+AP+CF" and "7+MAC+CF") codes. The results for the different 
codes are also labeled with 'A', 'B', and 'C, respectively. Since a MAC projection 
is not used in 7+AP+CF ('B'), no field energy is subtracted and it does not appear 
in the graph at lower-left. Note the relative size of increases in the kinetic energy of 
the codes, with smaller kinetic energies indicating stable oscillations. The eight- wave 
code with approximate projection and filter ("8+AP+CF") was run on this problem 
as well, but became unstable. 
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Fig. 11. Inclined flux tube results for 256, 512, and 1024 cells in each direction. 
These results used the 7- wave MHD code with filtering and a coefficient of 0.1. The 
quantities shown are the same as in Figure 9. 
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times smaller than the initial internal energy, and ~ 10 4 times smaller than 
the initial magnetic energy. Thus, small errors in either of these quantities can 
appear as sizeable changes in the kinetic energy. We note in particular that 
increases in the kinetic energy seem to correlate with increases in the fraction 
of the magnetic energy subtracted in either the filter or the MAC projection. 

The expected qualitative similarities between aligned and inclined flux tube 
results are present, such as oscillations in kinetic and magnetic energy of sim- 
ilar magnitudes. The size of the oscillations in the L 1 norm of v t and v n were 
also found to be very similar. However, substantial differences are also present. 
Comparing Figures 10 and 9, we note that the transient rise in kinetic energy 
early in the simulation begins later for the inclined tube than for the 1-D 
tube. In fact, the rise begins later and has smaller amplitude as resolution is 
increased, a fact that we attribute in part to the decrease in physical extent 
of the transition region. 

As a final test of the nonlinear behavior of the code, we reproduce the first 
MHD rotor problem outlined in TOO, earlier performed by Balsara and Spicer 
[5]. This problem constitutes a central, high-density (p = 10.0) region sur- 
rounded by a low-density medium. The central region is given a constant 
angular velocity. The entire domain is threaded by an initially constant field, 
B = 5.0e x . The rotation winds the magnetic field, sending Alfven waves prop- 
agating into the surrounding medium. 

This test was run using both seven-wave codes (MAC+CF and AP+CF) and 
one eight-wave variant (MAC+CF). We see no negative pressures in any of the 
code variants, and obtain first-order convergence (errors decreasing as the grid 
spacing to the 0.95 power, as averaged over all variables) using the seven-wave 
codes, and the eight-wave code with a MAC projection. A reproduction of 
Figure 18 from TOO is given in Figure 12. Note the steeper gradient in Mach 
number in comparison to the TOO result; we do not see the peaks therein 
attributed to pressure undershoots. 



5 Conclusions 

We have presented an unsplit method for ideal MHD which, when combined 
with projection and filtering steps, shows no effects of non-solenoidal fields, 
while retaining the co-location of all physical quantites at cell-centers. The 
latter point is important because such a uniform centering makes it easier to 
extend the scheme to adaptive meshes, and has the advantage of using a well- 
understood Godunov method for time integration. These are in contrast to the 
Constrained Transport (CT) approach, where staggered grids add additional 
software complexity. Furthermore, it is not obvious how to discretize diffusion 
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Fig. 12. Rotor problem results reproducing Figure 18 from TOO. The plot has the 
same number of contours and these contours lie between the same limits. A slightly 
expanded domain was used (x,y £ [—0.64,0.64]) in order to preserve the same grid 
spacing. 

operators to include non-ideal MHD effects on staggered grids. Lastly, we add 
that the effective advection scheme for the magnetic field components is not 
of the standard type, and it is unclear what its properties are in the presence 
of under-resolved gradients. 

The complexity and cost of the scheme relative to others is an important con- 
sideration. Our scheme requires three characteristic analysis steps plus twelve 
Riemann solution steps in 3-D. This is comparable to the six characteristic 
analysis and six Riemann solution steps required in a two-step Runge-Kutta 
scheme. The computational cost of projection is not insignificant, and adds 
additional software complexity. On a single grid and a single processor, the 
cost of the projection when computed using FFTs is a small fraction of the 
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overall cost of the computation (less than 200 floating-point operations per 
cell in 3D). Current research on analysis-based solvers for Poisson's equation 
[21,4] have the potential to make the cost of this remain less than that of the 
hyperbolic calculation, on both adaptive grids and on parallel processors. 

The Hodge projection of fields either at cell-edges (MAC projection) or cell- 
centers (in our case via an approximate projection) was found to be essential to 
accuracy and stability. The MAC projection was essential for accuracy of the 
eight-wave code in the presence of discontinuities. The seven-wave code results 
on the inclined shock tube were less sensitive to the type of projection used. 
The inclined flux tube results painted essentially the same picture, though the 
seven- wave MAC projected result was somewhat more robust. 

Use of a filter alone on linear problems reduces computational cost and does 
not affect accuracy. We find that both filter and projection are required on 
strongly nonlinear problems such as the inclined flux tube. Also essential to 
the accuracy of the scheme are the modifications suggested by Stone - linear 
waves not propagating along coordinate-axis directions are not second-order 
accurate without it. 

We found that in determining the accuracy of our code, both measuring the 
rate of convergence on nonlinear problems and comparing absolute magni- 
tude of the error were required. Significant differences were found between the 
convergence rates even of individual primitive variables. These data can be 
crucial in choosing between different algorithmic variants of a base code. We 
encourage future authors to incorporate convergence testing, taking care to 
calculate errors in each variable separately, of simple nonlinear problems such 
as the shock tubes run here into their test suites. 

Overall the ideal MHD code with a projection, either MAC or approximate, 
and a filter performed best on the suite of tests presented. Filtering of the 
magnetic field was important to code stability in the nonlinear problems, and 
accuracy in the deficient wave problems. 

The magnetized, perturbed flux tube constitutes a strong test of the stability 
of our schemes. The combination of a stationary discontinuity, low ratio of 
thermal to magnetic pressure (beta), and imposed perturbation caused sig- 
nificant problems for the base unsplit code without projection or filter. With 
a suitably chosen value of the filtering coefficient, the filter helped stability 
and worked to decrease the magnitude of a centered-difference measure of the 
divergence of the field. 
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A Appendix A: Divergence Constraints, Modified Equation Anal- 
ysis, and Eigenvector Deficiencies 



In this section, we will present a heuristic analysis of the effect of numeri- 
cal errors in the divergence-free constraint on the stability of finite-difference 
methods for the ideal MHD equations. Our starting point will be the modi- 
fied equation approach to analyzing the effect of truncation error on solution 
error. For any finite difference method, the modified equation is given by the 
original system of PDEs, with forcing terms given by the truncation error. In 
the present setting, the modified equation takes the following form: 

dtU Mod + y . pljjMod} = Tu (JJ Mod ) (A.l) 

V ■ B Mod = T D (U Mod ). (A.2) 

Here tjj is the usual truncation error for the numerical method obtained from 
applying the difference operator to a solution to the differential equation eval- 
uated on the grid. The equation for the evolution of To is obtained by taking 
the divergence of the modified equation for the evolution of U Mod , 

d t r D = V • t b . (A.3) 
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The truncation error forcing terms mimic the effect of numerical error on the 
computed solution. Specifically, we expect JJ Mod , the solution to the modified 
equation, to satisfy \\U Ax - U Mod \\ = 0(Ax p+1 ), where U Ax , the solution 
obtained from the p th -order scheme on a grid with spacing Ax, satisfies | \U Ax — 
U\\ = 0(Ax p ). 

In particular, for MHD, the effect of numerical error can be understood in 
terms of the truncation-error forcing in the modified equation causing the so- 
lution to violate the divergence-free constraint. Without that constraint being 
satisfied, the remaining ideal MHD equations can exhibit eigenvector defi- 
ciencies in the linearized-coeflicient matrix A, leading to anomalous loss of 
regularity and ill-posedness. In the numerical simulation, this translates into 
loss of accuracy and possibly instability of the underlying difference method. 

To see this, we consider the case of a small-amplitude wave corresponding to 
one of the eigenmodes of A (see Equation 8): 



W(x, t) = Wo + a(x - X k t)r k (A.4) 
A r fc = \ k r k r k = (f k , Of '. (A. 5) 

Then W(x, t) satisfies the MHD equations up to terms of 0(a 2 ). 

Without loss of generality, we take the direction of propagation to be in the 
x-direction in 3-D. However, we allow our computational spatial grid to have 
an arbitrary orientation in space. In that case, the modified equation corre- 
sponding to our numerical solution to the PDE in primitive variables W is 
given by (see Equation 7) 



dfW + A d x W Mod = r w (A.6) 
W Mod = (W Mod , B^ od ). (A. 7) 

If we define the new variable a Mod = l k ■ (W Mod — Wq), the modified equation 
dynamics can be reduced to the following system of two equations: 



d t a Mod + X k d x a Mod + (l k ■ a B )d x B^ od = l k ■ f (A.8) 

d t Bl lod = r B (A.9) 

In the system here, if one of the computational spatial grid coordinate axes 
is aligned with x, r B = 0. However, if the direction of propagation is not 
aligned with one of the computational coordinates, then in general t b ^ 0. 
If l k ■ a B 7^ and A = 0, then the left-hand side of A.8-A.9 is an example 
of a first-order system with an eigenvector deficiency. Such systems have an 
obvious loss of spatial regularity: a Mod grows like the derivative of r B . This 
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is in contrast to the behavior of well-posed hyperbolic systems, in which the 
solution has the same spatial regularity as the forcing. In other words, since 
regularity implies that there are as many derivatives in the solution as in the 
forcing term, discontinuous forcing of a hyperbolic system leaves the problem 
ill-posed. Discontinuity in tb implies that a Mod can grow without bound, owing 
to its dependence on the derivative of tb- 

In terms of a numerical method, we expect that the presence of such terms 
would lead to either an anomalous loss of accuracy or numerical instability. In 
the latter case, the forcing of a Mod in A. 8 takes the form of a finite difference 
operator applied to B^ Iod , whose spatial variation is due entirely to tb- If t~b 
fails to be smooth, either because of lack of smoothness in the initial data 
or in the finite difference formulae (eg. limiters), such lack of smoothness is 
immediately amplified. 

This discussion also provides an explanation for the behavior of the method 
described here. The use of the MAC projection and the filter does not elim- 
inate the truncation error terms that lead to the eigenvector deficiency, but 
regularizes it by smoothing. For example, the application of the filter in the 
plane-wave example corresponded to adding a diffusion term to the equation 
for Bf od , 

d t B™° d = T B + vd 2 x B™° d , (A. 10) 

with i] = O(Ax). The use of the filter alone is sufficient to stabilize the 
small-amplitude plane wave solution in Section 4.1, and in that case leads to 
a second-order accurate result. The MAC projection performs a more dras- 
tic smoothing, but only on the intermediate form of B^ od used to compute 
d x Bf od in equation A.8. 
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